/*
**  CXSC is a C++ library for eXtended Scientific Computing 
**
**  Copyright (C) 1990-2000 Institut fuer Angewandte Mathematik,
**                          Universitaet Karlsruhe, Germany
**            (C) 2000-2006 Wiss. Rechnen/Softwaretechnologie
**                          Universitaet Wuppertal, Germany
**
**  This library is free software; you can redistribute it and/or
**  modify it under the terms of the GNU Library General Public
**  License as published by the Free Software Foundation; either
**  version 2 of the License, or (at your option) any later version.
**
**  This library is distributed in the hope that it will be useful,
**  but WITHOUT ANY WARRANTY; without even the implied warranty of
**  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
**  Library General Public License for more details.
**
**  You should have received a copy of the GNU Library General Public
**  License along with this library; if not, write to the Free
**  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
*/

// MPI-Routines for C-XSC data types
// "Packed" version: Sends data in Packed format
//
// Author: Markus Grimmer
// Date: 5.3.2006
//
// Incorporated Data Types: real, interval, rvector, ivector, rmatrix, imatrix,
//                          complex, cinterval, cvector, civector, cmatrix, cimatrix,
//                          l_real, l_interval, l_complex,
//                          l_rvector, l_ivector, l_rmatrix, l_imatrix,
//                          dotprecision, idotprecision, 
//                          cdotprecision, cidotprecision
// Incorporated Communication Routines: MPI_Send, MPI_Isend, 
//                                      MPI_Recv

#include "cxsc_mpicomm.h"

using namespace cxsc;
using namespace std;

bool MPI_CXSC_TYPES_DEFINED = false;
MPI_Datatype MPI_CXSC_REAL;
MPI_Datatype MPI_CXSC_COMPLEX;
MPI_Datatype MPI_CXSC_INTERVAL;
MPI_Datatype MPI_CXSC_CINTERVAL;

// Definition of new MPI data types:

int MPI_Define_CXSC_Types(void){
	int err;
	if ((err = MPI_Type_contiguous(1,MPI_DOUBLE,&MPI_CXSC_REAL)) != MPI_SUCCESS){
		return err;
	}
	
	if ((err = MPI_Type_commit(&MPI_CXSC_REAL)) != MPI_SUCCESS){
		return err;
	}
	
	if ((err = MPI_Type_contiguous(2,MPI_CXSC_REAL,&MPI_CXSC_COMPLEX)) != MPI_SUCCESS){
		return err;
	}
	
	if ((err = MPI_Type_commit(&MPI_CXSC_COMPLEX)) != MPI_SUCCESS){
		return err;
	}
	
	if ((err = MPI_Type_contiguous(2,MPI_CXSC_REAL,&MPI_CXSC_INTERVAL)) != MPI_SUCCESS){
		return err;
	}
	
	if ((err = MPI_Type_commit(&MPI_CXSC_INTERVAL)) != MPI_SUCCESS){
		return err;
	}
	
	if ((err = MPI_Type_contiguous(2,MPI_CXSC_INTERVAL,&MPI_CXSC_CINTERVAL)) != MPI_SUCCESS){
		return err;
	}
	
	if ((err = MPI_Type_commit(&MPI_CXSC_CINTERVAL)) != MPI_SUCCESS){
		return err;
	}
	
	MPI_CXSC_TYPES_DEFINED = true;
	
	return err;
}

// =============================================================================
// MPI-Send/Recv for rmatrix:
// =============================================================================

int MPI_Send(rmatrix& rm, int dest, int tag, MPI_Comm MC){
	int err;
	
	if (!MPI_CXSC_TYPES_DEFINED){
		if ((err = MPI_Define_CXSC_Types()) != MPI_SUCCESS){
			return err;
		}
	}
	
	int bufferLength = sizeof(int) * 4, pos = 0;
	char *sendBuffer = new char[bufferLength];
	
	int lbRow = Lb(rm, ROW);
	int ubRow = Ub(rm, ROW);
	int lbCol = Lb(rm, COL);
	int ubCol = Ub(rm, COL);

	if((err = MPI_Pack(&lbRow, 1, MPI_INT, sendBuffer, bufferLength, &pos, MC)) != MPI_SUCCESS){
		return err;
	}
	if((err = MPI_Pack(&ubRow, 1, MPI_INT, sendBuffer, bufferLength, &pos, MC)) != MPI_SUCCESS){
		return err;
	}
	if((err = MPI_Pack(&lbCol, 1, MPI_INT, sendBuffer, bufferLength, &pos, MC)) != MPI_SUCCESS){
		return err;
	}
	if((err = MPI_Pack(&ubCol, 1, MPI_INT, sendBuffer, bufferLength, &pos, MC)) != MPI_SUCCESS){
		return err;
	}
	
	if((err = MPI_Send(sendBuffer, pos, MPI_PACKED, dest, tag, MC)) != MPI_SUCCESS){
		return err;
	}
	
	delete [] sendBuffer;
	sendBuffer = NULL;
	
	int numRows = ubRow - lbRow + 1;
	int numCols = ubCol - lbCol + 1;
	
	bufferLength = numRows * numCols * sizeof(double); 
	sendBuffer = new char[bufferLength];
	pos = 0;
	
	if((err = MPI_Pack(&rm[lbRow][lbCol], numRows * numCols, MPI_CXSC_REAL, sendBuffer, bufferLength, &pos, MC)) != MPI_SUCCESS){
		return err;
	}

	err = MPI_Send(sendBuffer, pos, MPI_PACKED, dest, tag, MC);
	
	delete [] sendBuffer;
	sendBuffer = NULL;
	
	return err;
}

int MPI_Recv(rmatrix& rm, int src, int tag, MPI_Comm MC, MPI_Status* MS){
	int err;
	
	if (!MPI_CXSC_TYPES_DEFINED){ 
		if ((err = MPI_Define_CXSC_Types()) != MPI_SUCCESS){
			return err;
		}
	}
	
	int pos = 0;
	int bufferLength = sizeof(int) * 4;
	char *recBuffer = new char[bufferLength];
	
	if ((err = MPI_Recv(recBuffer, bufferLength, MPI_PACKED, src, tag, MC, MS)) != MPI_SUCCESS){
		return err;
	}

	int lbRow = 0, ubRow = 0, lbCol = 0, ubCol = 0;
	
	if((err = MPI_Unpack(recBuffer, bufferLength, &pos, &lbRow, 1, MPI_INT, MC)) != MPI_SUCCESS){
		return err;
	}  
	if((err = MPI_Unpack(recBuffer, bufferLength, &pos, &ubRow, 1, MPI_INT, MC)) != MPI_SUCCESS){
		return err;
	}  
	if((err = MPI_Unpack(recBuffer, bufferLength, &pos, &lbCol, 1, MPI_INT, MC)) != MPI_SUCCESS){
		return err;
	}    
	if((err = MPI_Unpack(recBuffer, bufferLength, &pos, &ubCol, 1, MPI_INT, MC)) != MPI_SUCCESS){
		return err;
	}    
	
	delete [] recBuffer;
	recBuffer = NULL;
	
	int numRows = ubRow - lbRow + 1;
	int numCols = ubCol - lbCol + 1;
	
	Resize(rm, numRows, numCols);
	SetLb(rm, ROW, lbRow);
	SetLb(rm, COL, lbCol);
		
	bufferLength = numRows * numCols * sizeof(double);
	recBuffer = new char[bufferLength];

	pos = 0;
	
	if ((err = MPI_Recv(recBuffer, bufferLength, MPI_PACKED, src, tag, MC, MS)) != MPI_SUCCESS){
		return err;
	}
	
	if ((err = MPI_Unpack(recBuffer, bufferLength, &pos, &rm[lbRow][lbCol], numRows * numCols, MPI_CXSC_REAL, MC)) != MPI_SUCCESS){
		return err;
	}  
	
	delete [] recBuffer;
	recBuffer = NULL;
	
	return err;
}

int MPI_Bcast(rmatrix& rm, int rank, int root, MPI_Comm MC){
	int err;
	
	if (!MPI_CXSC_TYPES_DEFINED){
		if ((err = MPI_Define_CXSC_Types()) != MPI_SUCCESS){
			return err;
		}
	}
	
	int lbRow = 0, ubRow = 0, lbCol = 0, ubCol = 0;
	int bufferLength = sizeof(int) * 4, pos = 0;
	char *buffer = new char[bufferLength];
	
	if(rank == root){
		lbRow = Lb(rm, ROW);
		ubRow = Ub(rm, ROW);
		lbCol = Lb(rm, COL);
		ubCol = Ub(rm, COL);
		
		MPI_Pack(&lbRow, 1, MPI_INT, buffer, bufferLength, &pos, MC);
		MPI_Pack(&ubRow, 1, MPI_INT, buffer, bufferLength, &pos, MC);
		MPI_Pack(&lbCol, 1, MPI_INT, buffer, bufferLength, &pos, MC);
		MPI_Pack(&ubCol, 1, MPI_INT, buffer, bufferLength, &pos, MC);
	}
	
	if((err = MPI_Bcast(buffer, bufferLength, MPI_PACKED, root, MC)) != MPI_SUCCESS){
		return err;
	}
	
	if(rank != root){
		if(MPI_Unpack(buffer, bufferLength, &pos, &lbRow, 1, MPI_INT, MC) != MPI_SUCCESS){
			cout << "MPI_Unpack of lbRow failed." << endl;
		}  
		if(MPI_Unpack(buffer, bufferLength, &pos, &ubRow, 1, MPI_INT, MC) != MPI_SUCCESS){
			cout << "MPI_Unpack of ubRow failed." << endl;
		}  
		if(MPI_Unpack(buffer, bufferLength, &pos, &lbCol, 1, MPI_INT, MC) != MPI_SUCCESS){
			cout << "MPI_Unpack of lbCol failed." << endl;
		}    
		if(MPI_Unpack(buffer, bufferLength, &pos, &ubCol, 1, MPI_INT, MC) != MPI_SUCCESS){
			cout << "MPI_Unpack of ubCol failed." << endl;
		}
	}
	
	delete [] buffer;
	buffer = NULL;
		
	int numRows = ubRow - lbRow + 1;
	int numCols = ubCol - lbCol + 1;
			
	bufferLength = numRows * numCols * sizeof(double); 
	buffer = new char[bufferLength];
	pos = 0;
	
	if(rank == root){
		MPI_Pack(&rm[lbRow][lbCol], numRows * numCols, MPI_CXSC_REAL, buffer, bufferLength, &pos, MC);
	} else {
		Resize(rm, numRows, numCols);
		SetLb(rm, ROW, lbRow);
		SetLb(rm, COL, lbCol);
	}
	
	err = MPI_Bcast(buffer, bufferLength, MPI_PACKED, root, MC);
	
	if(rank != root){
		if(MPI_Unpack(buffer, bufferLength, &pos, &rm[lbRow][lbCol], numRows * numCols, MPI_CXSC_REAL, MC) != MPI_SUCCESS){
			cout << "MPI_Unpack of rm failed." << endl;
		}
	}
	
	delete [] buffer;
	buffer = NULL;
	
	return err;
}

int MPI_Scatter(rmatrix& m, rmatrix& sliceM, int sendcnt, int rank, int root, MPI_Comm MC){
	int err;
	
	if (!MPI_CXSC_TYPES_DEFINED){
		if ((err = MPI_Define_CXSC_Types()) != MPI_SUCCESS){
			return err;
		}
	}

	int *procBoundsSliceA = NULL, *boundsSliceA = NULL;

	if(rank == root){
		int slice = sendcnt;
		int commSize = 0;
		MPI_Comm_size(MPI_COMM_WORLD, &commSize);
		
		// Store the index bounds of matrix a.
		// Each index will be stored in the following order: [lbRow,ubRow,lbCol,ubCol, ...]
		boundsSliceA = new int[commSize * 4];

		for(int r = 0; r < commSize; r++){
			boundsSliceA[r*4]   = r * slice + 1;     // Rows lower bound
			boundsSliceA[r*4+1] = r * slice + slice; // Rows upper bound
			boundsSliceA[r*4+2] = Lb(m, COL);        // Cols lower bound
			boundsSliceA[r*4+3] = Ub(m, COL);        // Cols upper bound
		}
	}
	
	procBoundsSliceA = new int[4];
	
	MPI_Scatter(boundsSliceA, 4, MPI_INT, procBoundsSliceA, 4, MPI_INT, root, MPI_COMM_WORLD);

	delete [] boundsSliceA;
	boundsSliceA = NULL;
	
	int lbRow = procBoundsSliceA[0];
	int ubRow = procBoundsSliceA[1];
	int lbCol = procBoundsSliceA[2];
	int ubCol = procBoundsSliceA[3];

	delete [] procBoundsSliceA;
	procBoundsSliceA = NULL;
	
	char *sendBuffer = NULL;
	int sendBufferLength = 0;
	int numRows = 0, numCols = 0;
	int pos = 0;
	
	if(rank == root){
		numRows = Ub(m, ROW) - Lb(m, ROW) + 1;
		numCols = Ub(m, COL) - Lb(m, COL) + 1;
		
		sendBufferLength = numRows * numCols * sizeof(double);
		sendBuffer = new char[sendBufferLength];

		if((err = MPI_Pack(&m[lbRow][lbCol], numRows * numCols, MPI_CXSC_REAL, sendBuffer, sendBufferLength, &pos, MC)) != MPI_SUCCESS){
			return err;
		}
	}
	
	pos = 0;
	numRows = ubRow - lbRow + 1;
	numCols = ubCol - lbCol + 1;

	Resize(sliceM, numRows, numCols);
	SetLb(sliceM, ROW, lbRow);
	SetLb(sliceM, COL, lbCol);
	
	int recBufferLength = numRows * numCols * sizeof(double);
	char *recBuffer = new char[recBufferLength];

	err = MPI_Scatter(sendBuffer, recBufferLength, MPI_PACKED, recBuffer, recBufferLength, MPI_PACKED, root, MC); 
	
	if(MPI_Unpack(recBuffer, recBufferLength, &pos, &sliceM[lbRow][lbCol], numRows * numCols, MPI_CXSC_REAL, MC) != MPI_SUCCESS){
		cout << "MPI_Unpack of rm failed." << endl;
	}

	if(rank == root){
		delete [] sendBuffer;
		sendBuffer = NULL;
	}
	
	delete [] recBuffer;
	recBuffer = NULL;
	
	return err;
}

int MPI_Scatter(imatrix& m, imatrix& sliceM, int sendcnt, int rank, int root, MPI_Comm MC){
	int err;
	
	if (!MPI_CXSC_TYPES_DEFINED){
		if ((err = MPI_Define_CXSC_Types()) != MPI_SUCCESS){
			return err;
		}
	}

	int *procBoundsSliceA = NULL, *boundsSliceA = NULL;

	if(rank == root){
		int slice = sendcnt;
		int commSize = 0;
		MPI_Comm_size(MPI_COMM_WORLD, &commSize);
		
		// Store the index bounds of matrix a.
		// Each index will be stored in the following order: [lbRow,ubRow,lbCol,ubCol, ...]
		boundsSliceA = new int[commSize * 4];

		for(int r = 0; r < commSize; r++){
			boundsSliceA[r*4]   = r * slice + 1;     // Rows lower bound
			boundsSliceA[r*4+1] = r * slice + slice; // Rows upper bound
			boundsSliceA[r*4+2] = Lb(m, COL);        // Cols lower bound
			boundsSliceA[r*4+3] = Ub(m, COL);        // Cols upper bound
		}
	}
	
	procBoundsSliceA = new int[4];
	
	if((err = MPI_Scatter(boundsSliceA, 4, MPI_INT, procBoundsSliceA, 4, MPI_INT, root, MPI_COMM_WORLD)) != MPI_SUCCESS){
		return err;
	}

	delete [] boundsSliceA;
	boundsSliceA = NULL;
	
	int lbRow = procBoundsSliceA[0];
	int ubRow = procBoundsSliceA[1];
	int lbCol = procBoundsSliceA[2];
	int ubCol = procBoundsSliceA[3];

	delete [] procBoundsSliceA;
	procBoundsSliceA = NULL;
	
	char *sendBuffer = NULL;
	int sendBufferLength = 0;
	int numRows = 0, numCols = 0;
	int pos = 0;
	
	if(rank == root){
		numRows = Ub(m, ROW) - Lb(m, ROW) + 1;
		numCols = Ub(m, COL) - Lb(m, COL) + 1;
		
		sendBufferLength = numRows * numCols * sizeof(double) * 2;
		sendBuffer = new char[sendBufferLength];

		if((err = MPI_Pack(&m[lbRow][lbCol], numRows * numCols, MPI_CXSC_INTERVAL, sendBuffer, sendBufferLength, &pos, MC)) != MPI_SUCCESS){
			return err;
		}
	}
	
	pos = 0;
	numRows = ubRow - lbRow + 1;
	numCols = ubCol - lbCol + 1;

	Resize(sliceM, numRows, numCols);
	SetLb(sliceM, ROW, lbRow);
	SetLb(sliceM, COL, lbCol);
	
	int recBufferLength = numRows * numCols * sizeof(double);
	char *recBuffer = new char[recBufferLength];

	if((err = MPI_Scatter(sendBuffer, recBufferLength, MPI_PACKED, recBuffer, recBufferLength, MPI_PACKED, root, MC)) != MPI_SUCCESS) {
		return err;
	}
	
	if((err = MPI_Unpack(recBuffer, recBufferLength, &pos, &sliceM[lbRow][lbCol], numRows * numCols, MPI_CXSC_REAL, MC)) != MPI_SUCCESS){
		return err;
	}

	if(rank == root){
		delete [] sendBuffer;
		sendBuffer = NULL;
	}
	
	delete [] recBuffer;
	recBuffer = NULL;
	
	return err;
}

int MPI_Gather(rmatrix& m, rmatrix& sliceM, int root, MPI_Comm MC){
	int err;
	
	if (!MPI_CXSC_TYPES_DEFINED){
		if ((err = MPI_Define_CXSC_Types()) != MPI_SUCCESS){
			return err;
		}
	}
	
	int rank = 0;
	
	MPI_Comm_rank(MC, &rank);
	
	char *sendBuffer = NULL, *recBuffer = NULL;
	int sendBufferLength = 0, recBufferLength = 0;
	int numRows = 0, numCols = 0;
	int pos = 0;

	numRows = Ub(sliceM, ROW) - Lb(sliceM, ROW) + 1;
	numCols = Ub(sliceM, COL) - Lb(sliceM, COL) + 1;

	sendBufferLength = numRows * numCols * sizeof(double);
	sendBuffer = new char[sendBufferLength];

	if((err = MPI_Pack(&sliceM[Lb(sliceM, ROW)][Lb(sliceM, COL)], numRows * numCols, MPI_CXSC_REAL, sendBuffer, sendBufferLength, &pos, MC)) != MPI_SUCCESS){
		return err;
	}

	if(rank == root){
		numRows = Ub(m, ROW) - Lb(m, ROW) + 1;
		numCols = Ub(m, COL) - Lb(m, COL) + 1;
		
		recBufferLength = numRows * numCols * sizeof(double);
		recBuffer = new char[recBufferLength];
	}

	if((err = MPI_Gather(sendBuffer, sendBufferLength, MPI_PACKED, recBuffer, sendBufferLength, MPI_PACKED, root, MC)) != MPI_SUCCESS){
		return err;
	}

	delete [] sendBuffer;
	sendBuffer = NULL;
	
	if(rank == root){
		pos = 0;
		if((err = MPI_Unpack(recBuffer, recBufferLength, &pos, &m[Lb(m, ROW)][Lb(m, COL)], numRows * numCols, MPI_CXSC_REAL, MC)) != MPI_SUCCESS){
			return err;
		}

		delete [] recBuffer;
		recBuffer = NULL;
	}
		
	return err;
}

int MPI_Allgather(rmatrix& m, rmatrix& sliceM, MPI_Comm MC){
	int err;
	
	if (!MPI_CXSC_TYPES_DEFINED){
		if ((err = MPI_Define_CXSC_Types()) != MPI_SUCCESS){
			return err;
		}
	}

	char *sendBuffer = NULL, *recBuffer = NULL;
	int sendBufferLength = 0, recBufferLength = 0;
	int numRows = 0, numCols = 0;
	int pos = 0;

	numRows = Ub(sliceM, ROW) - Lb(sliceM, ROW) + 1;
	numCols = Ub(sliceM, COL) - Lb(sliceM, COL) + 1;

	sendBufferLength = numRows * numCols * sizeof(double);
	sendBuffer = new char[sendBufferLength];

	if((err = MPI_Pack(&sliceM[Lb(sliceM, ROW)][Lb(sliceM, COL)], numRows * numCols, MPI_CXSC_REAL, sendBuffer, sendBufferLength, &pos, MC)) != MPI_SUCCESS){
		return err;
	}

	numRows = Ub(m, ROW) - Lb(m, ROW) + 1;
	numCols = Ub(m, COL) - Lb(m, COL) + 1;
	
	recBufferLength = numRows * numCols * sizeof(double);
	recBuffer = new char[recBufferLength];

	if((err = MPI_Allgather(sendBuffer, sendBufferLength, MPI_PACKED, recBuffer, sendBufferLength, MPI_PACKED, MC)) != MPI_SUCCESS){
		return err;
	}

	delete [] sendBuffer;
	sendBuffer = NULL;
	
	pos = 0;
	if((err = MPI_Unpack(recBuffer, recBufferLength, &pos, &m[Lb(m, ROW)][Lb(m, COL)], numRows * numCols, MPI_CXSC_REAL, MC)) != MPI_SUCCESS){
		return err;
	}

	delete [] recBuffer;
	recBuffer = NULL;
	
	return err;
}

// =============================================================================
// MPI-Send/Recv for imatrix:
// =============================================================================

int MPI_Send(imatrix& im, int dest, int tag, MPI_Comm MC){
	int err;

	if (!MPI_CXSC_TYPES_DEFINED){ 
		if ((err = MPI_Define_CXSC_Types()) != MPI_SUCCESS){
			return err;
		}
	}

	int bufferLength = sizeof(int) * 4, pos = 0;
	char *sendBuffer = NULL;
	
	try {
		sendBuffer = new char[bufferLength];
	} catch(std::bad_alloc& e){
		cerr << e.what() << endl;
		return -1;
	}

	int lbRow = Lb(im, ROW);
	int lbCol = Lb(im, COL);
	int ubRow = Ub(im, ROW);
	int ubCol = Ub(im, COL);

	if((err = MPI_Pack(&lbRow, 1, MPI_INT, sendBuffer, bufferLength, &pos, MC)) != MPI_SUCCESS){
		return err;
	}
	if((err = MPI_Pack(&lbCol, 1, MPI_INT, sendBuffer, bufferLength, &pos, MC)) != MPI_SUCCESS){
		return err;
	}
	if((err = MPI_Pack(&ubRow, 1, MPI_INT, sendBuffer, bufferLength, &pos, MC)) != MPI_SUCCESS){
		return err;
	}
	if((err = MPI_Pack(&ubCol, 1, MPI_INT, sendBuffer, bufferLength, &pos, MC)) != MPI_SUCCESS){
		return err;
	}
	if((err = MPI_Send(sendBuffer, pos, MPI_PACKED, dest, tag, MC)) != MPI_SUCCESS){
		return err;
	}

	delete [] sendBuffer;
	sendBuffer = NULL;
	
	int numRows = ubRow - lbRow + 1;
	int numCols = ubCol - lbCol + 1;

	pos = 0;
	bufferLength = numRows * numCols * sizeof(double) * 2;

	try {
		sendBuffer = new char[bufferLength];
	} catch (std::bad_alloc& e){
		cerr << e.what() << endl;
		return -1;
	}

	if((err = MPI_Pack(&im[lbRow][lbCol], numRows * numCols, MPI_CXSC_INTERVAL, sendBuffer, bufferLength, &pos, MC)) != MPI_SUCCESS){
		return err;
	}

	err = MPI_Send(sendBuffer, pos, MPI_PACKED, dest, tag, MC);

	delete [] sendBuffer;
	sendBuffer = NULL;

	return err;
}

int MPI_Recv(imatrix& im, int src, int tag, MPI_Comm MC, MPI_Status* MS){
	int err;
	
	if (!MPI_CXSC_TYPES_DEFINED){ 
		if ((err = MPI_Define_CXSC_Types()) != MPI_SUCCESS){
			return err;
		}
	}
	
	int pos = 0;
	int bufferLength = sizeof(int) * 4;
	char *recBuffer = NULL;

	try {
		recBuffer = new char[bufferLength];
	} catch(std::bad_alloc& e){
		cerr << e.what() << endl;
		return -1;
	}

	if ((err = MPI_Recv(recBuffer, bufferLength, MPI_PACKED, src, tag, MC, MS)) != MPI_SUCCESS){
		return err;
	}

	int lbRow = 0, ubRow = 0, lbCol = 0, ubCol = 0;

	if ((err = MPI_Unpack(recBuffer, bufferLength, &pos, &lbRow, 1, MPI_INT, MC)) != MPI_SUCCESS){
		return err;
	}
	if ((err = MPI_Unpack(recBuffer, bufferLength, &pos, &lbCol, 1, MPI_INT, MC)) != MPI_SUCCESS){
		return err;
	}
	if ((err = MPI_Unpack(recBuffer, bufferLength, &pos, &ubRow, 1, MPI_INT, MC)) != MPI_SUCCESS){
		return err;
	}
	if ((err = MPI_Unpack(recBuffer, bufferLength, &pos, &ubCol, 1, MPI_INT, MC)) != MPI_SUCCESS){
		return err;
	}

	delete [] recBuffer;
	recBuffer = NULL;
	
	int numRows = ubRow - lbRow + 1;
	int numCols = ubCol - lbCol + 1;
	
	Resize(im, numRows, numCols);
	SetLb(im, ROW, lbRow);
	SetLb(im, COL, lbCol);

	pos = 0;
	bufferLength = numRows * numCols * sizeof(double) * 2;

	try {
		recBuffer = new char[bufferLength];
	} catch(std::bad_alloc& e){
		cerr << e.what() << endl;
		return -1;
	}

	if ((err = MPI_Recv(recBuffer, bufferLength, MPI_PACKED, src, tag, MC, MS)) != MPI_SUCCESS){
		return err;
	}

	err = MPI_Unpack(recBuffer, bufferLength, &pos, &im[lbRow][lbCol], numRows * numCols, MPI_CXSC_INTERVAL, MC);  

	delete [] recBuffer;
	recBuffer = NULL;

	return err;
}

int MPI_Gather(imatrix& m, imatrix& sliceM, int root, MPI_Comm MC){
	int err;
	
	if (!MPI_CXSC_TYPES_DEFINED){
		if ((err = MPI_Define_CXSC_Types()) != MPI_SUCCESS){
			return err;
		}
	}
	
	int rank = 0;
	
	MPI_Comm_rank(MC, &rank);

	char *sendBuffer = NULL, *recBuffer = NULL;
	int sendBufferLength = 0, recBufferLength = 0;
	int numRows = 0, numCols = 0;
	int pos = 0;

	numRows = Ub(sliceM, ROW) - Lb(sliceM, ROW) + 1;
	numCols = Ub(sliceM, COL) - Lb(sliceM, COL) + 1;

	sendBufferLength = numRows * numCols * sizeof(double) * 2;
	sendBuffer = new char[sendBufferLength];

	if((err = MPI_Pack(&sliceM[Lb(sliceM, ROW)][Lb(sliceM, COL)], numRows * numCols, MPI_CXSC_INTERVAL, sendBuffer, sendBufferLength, &pos, MC)) != MPI_SUCCESS){
		return err;
	}

	if(rank == root){
		numRows = Ub(m, ROW) - Lb(m, ROW) + 1;
		numCols = Ub(m, COL) - Lb(m, COL) + 1;
		
		recBufferLength = numRows * numCols * sizeof(double) * 2;
		recBuffer = new char[recBufferLength];
	}

	if((err = MPI_Gather(sendBuffer, sendBufferLength, MPI_PACKED, recBuffer, sendBufferLength, MPI_PACKED, root, MC)) != MPI_SUCCESS){
		return err;
	}

	delete [] sendBuffer;
	sendBuffer = NULL;
	
	if(rank == root){
		pos = 0;
		if((err = MPI_Unpack(recBuffer, recBufferLength, &pos, &m[Lb(m, ROW)][Lb(m, COL)], numRows * numCols, MPI_CXSC_INTERVAL, MC)) != MPI_SUCCESS){
			return err;
		}

		delete [] recBuffer;
		recBuffer = NULL;
	}
		
	return err;
}
